MODULE smmboot
  USE GLOBVAR
  USE NRUTIL
  USE PROBABILITY
  USE MATRIX
  USE STATISTICS
  USE random
  IMPLICIT NONE
  
    CONTAINS
   REAL(8) FUNCTION momb(thetain) 
   !SUBROUTINE mom(theta_ini,beta_f,check,N_int,sc)
        IMPLICIT NONE
        REAL(8), INTENT(IN):: thetain(:) !initial guess

        REAL(8):: Ls(N_inds,N_testc),Ep(N_inds,N_testc),as,pen,pen2
        REAL(8):: L0s(N_inds),eta(N_inds,N_testc),e_L0(n_inds), e_eta(n_inds),mom_simu(N_mom),dif(N_mom,1)!, latent(N_inds*N_testl,N_test),ps(N_testl,N_test)
        REAL(8):: inters(n_inds,n_levelc,20,3),intercs(n_inds,n_levelc,3)
        REAL(8),ALLOCATABLE:: AVSE(:),cors(:,:)
        INTEGER:: ducs(n_inds,n_levelc),dus(n_levelc),inds(N_inds,n_levelc) 
        ! parameters
        REAL(8):: beta_mu(6,1),beta_eta(n_x_eta,1),var_l,Cbar(N_testc),var_e(N_testc), etav,delta(N_levelc),beta_y(2,1)!,beta_level(n_l,1)
        INTEGER:: Ds(N_inds,N_testc),indexc(n_testc,2), artcs(n_inds,n_levelc,11) 
        INTEGER:: i,j,k,jj,kk,m,t,n,l,tt,mmm
        
        !!!!!!!!!!!!!!!!!!Simulate Moments
        CALL SET_SEED(1438914432,727629,19890414,111111)
        beta_mu(1:n_x_l0,1)=thetain(1:n_x_l0)                                              ! Initial Latent skill
        DO i=2,4
          beta_mu(i,1)=1.40d0**(thetain(i))/(10.0d0*(1.0d0+1.40d0**(thetain(i))))
        END DO
        beta_mu(5,1)=1.0d0-1.40d0**(thetain(5))
        
        var_l=4.0d0*(1.4d0**thetain(n_x_l0+1)/(1.0d0+1.4d0**thetain(n_x_l0+1)))     ! variance of the shoch of initial latent skill
                Cbar=1.0d0
        Cbar(2:N_levelc) =thetain(n_x_l0+2:n_x_l0+N_levelc)                ! mimimum requriement for each difficulty level
        var_e=1.0d0
        var_e(2:N_levelc-1)=1.4d0**thetain(1+n_x_l0+N_levelc:n_x_l0+2*N_levelc-2)     ! variance of the shock at each difficulty level
        beta_eta(1:n_x_eta,1)=thetain(n_x_l0+2*N_levelc-1:n_x_l0+2*N_levelc+n_x_eta-2)           ! x for eta
        etav=1.40d0**thetain(n_x_l0+2*N_levelc+n_x_eta-1)                           ! variance of shock for eta
        delta=1.0d0
        delta(2)=4*(1.10d0**thetain(n_x_l0+2*N_levelc+n_x_eta))/(1.0d0+1.10d0**thetain(n_x_l0+2*N_levelc+n_x_eta))  ! delta for each level
        DO i=5,8
            delta(i)=4*(1.10d0**thetain(n_x_l0+2*N_levelc+n_x_eta+i-4))/(1.0d0+1.10d0**thetain(n_x_l0+2*N_levelc+n_x_eta+i-4))  ! delta for each level
        END DO
        DO i=10,n_levelc-1
            delta(i)=4*(1.10d0**thetain(n_x_l0+2*N_levelc+n_x_eta+i-5))/(1.0d0+1.10d0**thetain(n_x_l0+2*N_levelc+n_x_eta+i-5))  ! delta for each level
        END DO
                beta_y=0.0d0
                     beta_y(1:2,1)=thetain(n_x_l0+3*N_levelc+n_x_eta-6+1:n_x_l0+3*N_levelc+n_x_eta-6+2)
        !        beta_level=0.0d0
        !beta_level(2:n_l,1)=thetain(n_x_l0+3*N_levelc+n_x_eta+2:n_x_l0+3*N_levelc+n_x_eta+n_l)

        ! generate shock terms
        
        DO i=1,N_inds
            e_L0(i) =Sample_NORMAL(0.0d0,var_l)
            k=i-INT((i-1)/N_IND)*N_ind
            e_eta(i)=Sample_Normal(intermean(mmb(k)),etav)
        END DO
        
        ! GENERATE INITIAL LATENT VALUES AND 
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            L0s(i)=1.40d0**(DOT_PRODUCT(L0(mmb(k),:),beta_mu(:,1))+e_L0(i))
        END DO
        
       ! SIMULATE LATENT SKILL PROCESS ! HERE ONLY FOR LANGUAGE SKILLS 
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testc
                IF (j<firstc(mmb(k),1)) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j==firstc(mmb(k),1) .and. m_c(mmb(k),j)==1) THEN
                    IF (x_testc(mmb(k),j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testc(mmb(k),j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testc(mmb(k),j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    
                    kk=levelc(j)
68                  Ls(i,j)=L0s(i)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthc(j)+beta_y(2,1)*DBLE(timec(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 68 
                ELSE IF (j==firstc(mmb(k),1) .and. m_c(mmb(k),j)==-1) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j>firstc(mmb(k),1) .and. Ls(i,j-1)<0.0d0 .and. m_c(mmb(k),j)==-1) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j>firstc(mmb(k),1) .and. Ls(i,j-1)<0.0d0 .and. m_c(mmb(k),j)==1) THEN
                    IF (x_testc(mmb(k),j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testc(mmb(k),j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testc(mmb(k),j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levelc(j)
77                  Ls(i,j)=L0s(i)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthc(j)+beta_y(2,1)*DBLE(timec(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 77
                ELSE IF (j>firstc(mmb(k),1) .and. Ls(i,j-1)>=0.0d0 .and. m_c(mmb(k),j)==-1 ) THEN
                    Ls(i,j)=Ls(i,j-1)
                ELSE IF (j>firstc(mmb(k),1) .and. Ls(i,j-1)>=0.0d0 .and. m_c(mmb(k),j)==1) THEN
                    IF (x_testc(mmb(k),j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testc(mmb(k),j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testc(mmb(k),j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levelc(j)
84                  Ls(i,j)=Ls(i,j-1)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthc(j)+beta_y(2,1)*DBLE(timec(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 84 
                END IF
            END DO
        END DO
        
        ! GENERATE TASK PERFORMANCE !
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testc
                kk=levelc(j)
                IF (m_c(mmb(k),j)==-1) THEN
                    DS(i,j)=-99
                ELSE IF (Ls(i,j)>=SUM(Cbar(1:kk)) .and. m_c(mmb(k),j)==1) THEN
                    DS(i,j)=1
                ELSE IF (Ls(i,j)<SUM(Cbar(1:kk)) .and. m_c(mmb(k),j)==1) THEN
                    DS(i,j)=0
                END IF
            END DO
        END DO
                
         !!!!!!!!!!!!!!!!!From here !!!!!!!!!!!!!!!!!!!!!!!
        
    !!! cognitive tasks!!!!
    indexc=0
    DO i=1,N_inds
      DO j=1,N_testc
         IF (ds(i,j)>-90) THEN
             indexc(j,1)=indexc(j,1)+1
             IF (ds(i,j)>0.5) indexc(j,2)=indexc(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testc
        IF (indexc(i,1)>0.5) THEN
            mom_simu(i)=DBLE(indexc(i,2))/DBLE(indexc(i,1)) ! PASSING RATE FOR EACH Cognitive TASK
        ELSE 
            mom_simu(i)=-99
        END IF
    END DO
    
    !!!! generate artc variables !!!!
    artcs=-1
    DO i=1,N_inds
        DO j=1,N_testc
            IF (ds(i,j)>-90) THEN
                kk=levelc(j)
                k=i-INT((i-1)/N_IND)*N_ind
                m=repc(mmb(k),j)
                artcs(i,kk,m)=ds(i,j)
                inters(i,kk,m,1:3)=x_testc(mmb(k),j,1:3)
            END IF
        END DO
    END DO
    
    !!!! passing rate at each level!!!!

    DO j=N_testc+1,N_testc+N_levelc
       m=0
       n=0 
       DO i=1,N_inds
           DO k=1,11
               IF (artcs(i,j-N_testc,k)>-0.5) THEN
                   m=m+1
                   IF (artcs(i,j-N_testc,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           mom_simu(j)=DBLE(n)/DBLE(m)
       ELSE
           mom_simu(j)=-99
       END IF
    END DO
    
    !!!! correlation across levels !!!!
    
    mmm=N_testc+N_levelc
    
    l=0
    DO j=1,N_levelc-2
        IF (j<=2 .or. (j>4 .and. j<=8)) THEN
            tt=5
        ELSE IF (j>2 .and. j<5) THEN
            tt=2
        ELSE IF (j==9) THEN
            tt=2
        ELSE IF (j==10) THEN
            tt=4
        ELSE IF (j==11) THEN
            tt=5
        END IF
        
        IF (j<2 .or. (j>3 .and. j<8)) THEN
            kk=5
        ELSE IF (j>1 .and. j<4) THEN
            kk=2
        ELSE IF (j==8) THEN
            kk=2
        ELSE IF (j==9) THEN
            kk=4
        ELSE IF (j==10) THEN
            kk=5
        ELSE IF (j==11) THEN
            kk=3
        END IF
        
         DO t=1,tt
           DO k=1,kk
             m=0
             n=0 
             l=l+1
             DO i=1,N_inds 
                IF (artcs(i,j,t)==1 .and. artcs(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                    m=m+1
                    IF (artcs(i,j,t)==1 .and. artcs(i,j+1,k)==1) n=n+1
                END IF
             END DO
             
             IF (m>0) THEN
                mom_simu(mmm+L)=DBLE(n)/DBLE(m)
             ELSE
                mom_simu(mmm+L)=-99
             END IF
           END DO
         END DO
    END DO
    
    mmm=N_testc+N_levelc+177
    
    l=0
    DO j=1,N_levelc-1
        IF (j<=2 .or. (j>4 .and. j<=8)) THEN
            tt=4
            kk=5
        ELSE IF (j>2 .and. j<5) THEN
            tt=1
            kk=2
        ELSE IF (j==9) THEN
            tt=1
            kk=2
        ELSE IF (j==10) THEN
            tt=3
            kk=4
        ELSE IF (j==11) THEN
            tt=4
            kk=5
        ELSE IF (j==12) THEN
            tt=2
            kk=3
        END IF
        
         DO t=1,tt
           DO k=t+1,kk
             m=0
             n=0 
             l=l+1
             DO i=1,N_inds 
                IF (artcs(i,j,t)==1 .and. artcs(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                    m=m+1
                    IF (artcs(i,j,t)==1 .and. artcs(i,j,k)==1) n=n+1
                END IF
             END DO
             
             IF (m>0) THEN
                mom_simu(mmm+L)=DBLE(n)/DBLE(m)
             ELSE
                mom_simu(mmm+L)=-99
             END IF
           END DO
         END DO
    END DO
    
     mmm=N_testc+N_levelc+177+82
    
    l=0
    DO j=1,N_levelc-3
        IF (j==1) THEN
            tt=5
            kk=2
        ELSE IF (j==2) THEN
            tt=5
            kk=2
        ELSE IF (j==3) THEN
            tt=2
            kk=5
        ELSE IF (j==4) THEN
            tt=2
            kk=5
        ELSE IF (j==5) THEN
            tt=5
            kk=5
        ELSE IF (j==6) THEN
            tt=5
            kk=5
        ELSE IF (j==7) THEN
            tt=5
            kk=2
        ELSE IF (j==8) THEN
            tt=5
            kk=4
        ELSE IF (j==9) THEN
            tt=2
            kk=5
        ELSE IF (j==10) THEN
            tt=4
            kk=3
        END IF
        
        
      DO t=1,tt
        DO k=1,kk
          m=0
          n=0 
          l=l+1
          DO i=1,N_inds 
             IF (artcs(i,j,t)==1 .and. artcs(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artcs(i,j,t)==1 .and. artcs(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
        END DO
      END DO
    END DO  
    
     mmm=N_testc+N_levelc+177+82+142
    
    
    indexc=0
    DO i=1,N_inds
         k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testc
          as=firstc(mmb(k),2)-monthc(j)-DBLE(timec(j))/4.0d0
         IF (ds(i,j)>-90 .and. as>-1.25d0) THEN
             indexc(j,1)=indexc(j,1)+1
             IF (ds(i,j)>0.5) indexc(j,2)=indexc(j,2)+1
         END IF
      END DO
    END DO
    
    DO i=1,N_testc
        IF (indexc(i,1)>0.5) THEN
            mom_simu(mmm+i)=DBLE(indexc(i,2))/DBLE(indexc(i,1)) ! PASSING RATE FOR EACH cognitive TASK
        ELSE 
            mom_simu(mmm+i)=-99
        END IF
    END DO
    
    mmm=2*N_testc+N_levelc+177+82+142
    
    indexc=0
    DO i=1,N_inds
         k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testc
          as=DBLE(firstc(mmb(k),2)-monthc(j)-DBLE(timec(j))/4.0d0)
         IF (Ds(i,j)>-90 .and. as<-1.25d0) THEN
             indexc(j,1)=indexc(j,1)+1
             IF (ds(i,j)>0.5) indexc(j,2)=indexc(j,2)+1
         END IF
      END DO
    END DO
    
    DO i=1,N_testc
        IF (indexc(i,1)>0.5) THEN
            mom_simu(mmm+i)=DBLE(indexc(i,2))/DBLE(indexc(i,1)) ! PASSING RATE FOR EACH cognitive TASK
        ELSE 
            mom_simu(mmm+i)=-99
        END IF
    END DO
    
    mmm=3*N_testc+N_levelc+177+82+142
    
    l=0
    DO j=1,N_levelc-1
      DO t=1,5
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artcs(i,j,t)>-0.5) THEN   
                 m=m+1
                 IF (artcs(i,j,t)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(mmm+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(mmm+l)=-99
          END IF
      END DO
    END DO
    mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)  
    
    DUcs=-99
    inds=-99
    DO i=1,N_inds
        DO k=1,N_levelc
            DO j=1,11
                IF (artcs(i,k,j)==1 .and. j==1) THEN
                    ducs(i,k)=1
                    inds(i,k)=1
                ELSE IF (artcs(i,k,j)==0 .and. j==1) THEN
                    ducs(i,k)=j
                    inds(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF( artcs(i,k,j-1)==0 .and. artcs(i,k,j)==0 .and. inds(i,k)==0) THEN
                       ducs(i,k)=j
                       inds(i,k)=0
                       IF (j<11 ) THEN
                           IF (artcs(i,k,j+1)==-1) THEN
                           ducs(i,k)=j+1
                           inds(i,k)=1
                           END IF
                       END IF 
                    ELSE IF (j>1 .and. artcs(i,k,j-1)==0 .and. artcs(i,k,j)==0 .and. inds(i,k)==1) THEN
                        ducs(i,k)=ducs(i,k)
                        inds(i,k)=1
                    ELSE IF (j>1 .and. artcs(i,k,j-1)==0 .and. artcs(i,k,j)==1 .and. inds(i,k)==0) THEN
                        ducs(i,k)=j
                        inds(i,k)=1
                    ELSE IF (j>1 .and. artcs(i,k,j-1)==0 .and. artcs(i,k,j)==1 .and. inds(i,k)==1) THEN
                        ducs(i,k)=ducs(i,k)
                        inds(i,k)=1   
                    END IF
                END IF
                
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levelc
        m=0
        n=0
        DO i=1,n_inds
            IF (ducs(i,k)>-90) THEN
               m=m+1
               n=n+ducs(i,k)
            END IF
        END DO
        IF (m>0) THEN
          mom_simu(mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           mom_simu(mmm+k)=-99
        END IF
    END DO
    
    mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)+n_levelc 
    
        dus=0
    DO k=1,n_levelc
        DO i=1,N_inds
            IF (ducs(i,k)>-90) THEN
                dus(k)=dus(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levelc
        allocate(avse(dus(k)))
        m=0
        DO i=1,n_inds
            IF (ducs(i,k)>-90) THEN
                m=m+1
                avse(m)=ducs(i,k)
            END IF
        END DO
        mom_simu(mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
      mmm=3*N_testc+N_levelc+177+82+142+5*(N_levelc-1)+2*n_levelc 

    INDs=0
    intercs=0.0d0
    DO k=1,N_levelc
        DO i=1,N_inds
          DO j=1,11
              IF (inters(i,k,j,1)>-90) THEN
                  inds(i,k)=inds(i,k)+1
                  intercs(i,k,1)=intercs(i,k,1)+inters(i,k,j,1)
                  intercs(i,k,2)=intercs(i,k,2)+inters(i,k,j,2)
                  intercs(i,k,3)=intercs(i,k,3)+inters(i,k,j,3)
              END IF
          END DO
          IF (inds(i,k)>0) THEN
               intercs(i,k,1)=DBLE(intercs(i,k,1))/DBLE(inds(i,k))
               intercs(i,k,2)=DBLE(intercs(i,k,2))/DBLE(inds(i,k))
               intercs(i,k,3)=DBLE(intercs(i,k,3))/DBLE(inds(i,k))
          END IF
        END DO
    END DO
    
    dus=0
    DO k=1,n_levelc
        DO i=1,N_inds
            IF (inds(i,k)>0 .AND. ducs(i,k)>-90) THEN
            dus(k)=dus(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levelc
        allocate(cors(dus(k),4))
        m=0
        DO i=1,n_inds
             IF (inds(i,k)>0 .AND. ducs(i,k)>-90) THEN
                m=m+1
                cors(m,1)=intercs(i,k,1)
                cors(m,2)=intercs(i,k,2)
                cors(m,3)=intercs(i,k,3)
                cors(m,4)=ducs(i,k)
             END IF
        END DO
        IF (k<n_levelc) THEN
        mom_simu(mmm+3*(k-1)+1)=Correlation(cors(:,1),cors(:,4))
        mom_simu(mmm+3*(k-1)+2)=Correlation(cors(:,2),cors(:,4))
        mom_simu(mmm+3*(k-1)+3)=Correlation(cors(:,3),cors(:,4))
        ELSE 
            mom_simu(mmm+3*(k-1)+1)=-99
            mom_simu(mmm+3*(k-1)+2)=-99
            mom_simu(mmm+3*(k-1)+3)=-99
        END IF
        
            
        DEALLOCATE(cors)
    END DO

    
    
    !    
    !     !!!!!!!!!!!!! Now finish simulate the fake data and start to generat the moments !!!!!!!!!!!!!!!!!!!!!!!!!
    !     
    
    DO i=1,N_mom
        IF (ABS(moment(i))<1.0d-10) THEN
            dif(i,1)=((mom_simu(i)-moment(i)))**2.0D0
        ELSE 
            dif(i,1)=((DBLE((mom_simu(i)-moment(i))))**2.0D0/DBLE(ABS(moment(i))))
        END IF 
        IF (i>N_testc+N_levelc) THEN
           dif(i,1)=1.0d0*dif(i,1)
        END IF 
        !IF (mom_simu(i)>-10.0d0) THEN        
        !      pen=exp(5.0d0*((1.0d0/(1.009d0-mom_simu(i)))/100.0d0))
        !      dif(i,1)=pen*dif(i,1)
        !END IF
        !IF (mom_simu(i)>-10.0d0) THEN        
        !      pen2=exp(5.0d0*((1.0d0/(0.009d0+mom_simu(i)))/100.0d0))
        !      dif(i,1)=pen2*dif(i,1)
        !END IF

      
        !dif(i)=(100.0d0/LOG(ABS(dif(i))))*dif(i)
        !dif(i)=dif(i)**2.0d0
        !IF (i<=25) dif(i)=dif(i)*25.0d0
    END DO
    
    !mm=MATMUL(MATMUL(TRANSPOSE(DIF),WEI),DIF)
    momb=SUM(dif)
    
    
    
                OPEN(6424,file="dif.out")
                  DO i = 1,N_mom
                    WRITE(6424,'(200F16.8)') SQRT(dif(i,1)),mom_simu(i),moment(i)
                  END DO
                CLOSE(6424)
            
   
   END FUNCTION momb
    
       
    
END MODULE smmboot